import logging
import os
from pathlib import Path
from IPython.display import display
from PIL import Image

# 调整参数
denose_nb_neighbors = 80  # 噪声的相邻要素数
denose_std_ratio = 0.5     # 噪声去除的标准偏差比率
# resolution = 0.1           # DEM分辨率，值越小越精细但慢（例如 0.5=每个像素显示0.5平方米的高程值）
# threshold = 0.7            # 高程变化阈值（以米为单位）


def get_las_files(base_dir, date_str):
    """获取指定日期目录下的所有.las文件"""
    logging.info(f"获取 {date_str} 目录下的所有.las文件")
    dir_path = Path(base_dir) / f'{date_str}'
    if not dir_path.exists():
        raise FileNotFoundError(f"目录不存在: {dir_path}")
    return sorted(str(path) for path in dir_path.glob('*.las'))




# 剔除离群孤立点
def denoise_las(input_path, nb_neighbors, std_ratio):
    import laspy
    import open3d as o3d
    import numpy as np

    las = laspy.read(input_path)
    points = np.vstack((las.x, las.y, las.z)).transpose()
    logging.info(f'已加载 {input_path}')

    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(points)

    # 剔除孤立点
    cl, ind = pcd.remove_statistical_outlier(nb_neighbors=nb_neighbors,
                                             std_ratio=std_ratio)
    filtered = pcd.select_by_index(ind)

    logging.info(f'已剔除孤立点 {input_path}')

    return np.asarray(filtered.points)



import numpy as np
import rasterio
from rasterio.transform import from_origin
from scipy.interpolate import griddata

def generate_dem_2(point_arrays, res, dem_path=None):

    pts = np.vstack(point_arrays)
    x, y, z = pts[:,0], pts[:,1], pts[:,2]
    x_min_in, x_max_in = x.min(), x.max()
    y_min_in, y_max_in = y.min(), y.max()
    logging.info(f"输入点云 X 范围: [{x_min_in}, {x_max_in}], Y 范围: [{y_min_in}, {y_max_in}]")

    if x_min_in < 0 or y_min_in < 0:
        logging.info("警告：输入点云存在负坐标，可能是投影不当或坐标轴混淆。请检查 point_arrays。")

    xmin, xmax = x_min_in, x_max_in
    ymin, ymax = y_min_in, y_max_in
    if xmin >= xmax or ymin >= ymax:
        raise ValueError("点云坐标范围异常：xmin>=xmax 或 ymin>=ymax，请检查原始数据。")

    nx = int(np.ceil((xmax - xmin) / res))
    ny = int(np.ceil((ymax - ymin) / res))

    xi = np.linspace(xmin, xmax, nx)
    yi = np.linspace(ymin, ymax, ny)
    grid_x, grid_y = np.meshgrid(xi, yi)

    if not np.isclose(grid_x.min(), x_min_in) or not np.isclose(grid_x.max(), x_max_in):
        logging.info(f"grid_x 范围 [{grid_x.min():.6f}, {grid_x.max():.6f}] 与输入 [{x_min_in:.6f}, {x_max_in:.6f}] 不一致")
    if not np.isclose(grid_y.min(), y_min_in) or not np.isclose(grid_y.max(), y_max_in):
        logging.info(f"grid_y 范围 [{grid_y.min():.6f}, {grid_y.max():.6f}] 与输入 [{y_min_in:.6f}, {y_max_in:.6f}] 不一致")

    grid_z = griddata((x, y), z, (grid_x, grid_y), method='nearest')

    if dem_path:
        logging.info(f"写入 GeoTIFF 文件：{dem_path}")
        transform = from_origin(xmin, ymax, res, res)
        with rasterio.open(
            dem_path, 'w', driver='GTiff',
            height=ny, width=nx, count=1,
            dtype=grid_z.dtype,
            crs='+proj=utm +zone=48 +datum=WGS84',
            transform=transform
        ) as dst:
            dst.write(grid_z, 1)

    return grid_x, grid_y, grid_z




# 对两个时间点的坐标交集，计算出高程差（第一时间点减第二时间点）
import numpy as np
from scipy.interpolate import griddata

def dem_diff_overlap(gx1, gy1, gz1, gx2, gy2, gz2):
    xmin = max(gx1.min(), gx2.min())
    xmax = min(gx1.max(), gx2.max())
    ymin = max(gy1.min(), gy2.min())
    ymax = min(gy1.max(), gy2.max())

    if xmin >= xmax or ymin >= ymax:
        raise ValueError("两 DEM 无重叠区域")

    res_x = gx1[0,1] - gx1[0,0]
    res_y = gy1[1,0] - gy1[0,0]

    xi = np.arange(xmin, xmax + res_x/2, res_x)
    yi = np.arange(ymin, ymax + res_y/2, res_y)
    grid_x3, grid_y3 = np.meshgrid(xi, yi)

    pts1 = np.vstack([gx1.ravel(), gy1.ravel()]).T
    vals1 = gz1.ravel()
    gz1_res = griddata(pts1, vals1, (grid_x3, grid_y3), method='nearest')

    pts2 = np.vstack([gx2.ravel(), gy2.ravel()]).T
    vals2 = gz2.ravel()
    gz2_res = griddata(pts2, vals2, (grid_x3, grid_y3), method='nearest')

    gz_diff = gz1_res - gz2_res

    return grid_x3, grid_y3, gz_diff




def filter_by_realland_mask(gx, gy, gz, mask_path):

    with rasterio.open(mask_path) as msrc:
        bands = msrc.read()
        tr = msrc.transform
        Hm, Wm = bands.shape[1], bands.shape[2]

    alpha = bands[3] if bands.shape[0] == 4 else np.full((Hm, Wm), 255, dtype=np.uint8)

    xs = gx.ravel()
    ys = gy.ravel()
    dz = gz.ravel()
    valid = ~np.isnan(dz)
    xs, ys, dz = xs[valid], ys[valid], dz[valid]

    a, b, c, d, e, f = tr.a, tr.b, tr.c, tr.d, tr.e, tr.f
    cols = np.floor((xs - c) / a + 0.5).astype(int)
    rows = np.floor((ys - f) / e + 0.5).astype(int)

    in_bounds = (rows>=0) & (rows < Hm) & (cols>=0) & (cols < Wm)
    rows, cols, xs, ys, dz = rows[in_bounds], cols[in_bounds], xs[in_bounds], ys[in_bounds], dz[in_bounds]

    mask = alpha[rows, cols] != 0
    rows, cols, xs, ys, dz = rows[mask], cols[mask], xs[mask], ys[mask], dz[mask]

    return rows, cols, xs, ys, dz

def overlay_blue_on_rgb(xs, ys, rgb_path, out_path=None):

    with rasterio.open(rgb_path) as src:
        arr = src.read()
        tr  = src.transform
        Hb, Wb = src.height, src.width
        logging.info(tr)

    # 反算到 rgb 行列
    a,b,c,d,e,f = tr.a, tr.b, tr.c, tr.d, tr.e, tr.f
    # 计算右下角的地理坐标
    x_right_bottom = a * (Wb - 1) + b * (Hb - 1) + c
    y_right_bottom = d * (Wb - 1) + e * (Hb - 1) + f

    cols = np.floor((xs - c)/a + 0.5).astype(int)
    rows = np.floor((ys - f)/e + 0.5).astype(int)


    # 读底图并转 PIL
    img = reshape_as_image(arr)
    pil = Image.fromarray(img).convert("RGBA")
    overlay = Image.new("RGBA", (Wb, Hb), (0,0,0,0))
    draw = ImageDraw.Draw(overlay)
    for r, c in zip(rows, cols):
        draw.rectangle([(c, r),(c+1, r+1)], fill=(0,0,255,255))

    comp = Image.alpha_composite(pil, overlay)

    if out_path:
        logging.info(f"写入标记下塌图片：{out_path}")
        comp.save(out_path)

    # 返回标记后的图像，左上角和右下角的地理坐标
    return comp, tr.c, tr.f, x_right_bottom, y_right_bottom






import numpy as np
import rasterio
from rasterio.plot import reshape_as_image
from rasterio.transform import xy
from PIL import Image, ImageDraw, ImageFont
import cv2

def mark_color_blocks_with_number_2(pil_image,
                                   target_rgb=(0,0,255), tol=10,
                                   min_area=500, max_area=10000,
                                   outline_color=(255,0,0), outline_width=2,
                                   text_color=(255,0,0), font_size=20,
                                   center_dist_thresh=5, font_path="DejaVuSans.ttf"):
    W, H = pil_image.size
    img = np.array(pil_image.convert("RGB"))
    lower = np.clip(np.array(target_rgb)-tol,0,255).astype(np.uint8)
    upper = np.clip(np.array(target_rgb)+tol,0,255).astype(np.uint8)
    mask = cv2.inRange(img, lower, upper)
    num, labels, stats, _ = cv2.connectedComponentsWithStats(mask,8)
    draw = ImageDraw.Draw(pil_image)
    try:
        font = ImageFont.truetype(font_path, font_size)
    except:
        font = ImageFont.load_default()
    centers = []; center_pix = []; idx=1
    for i in range(1, num):
        x,y,w,h,area = stats[i]
        if area<min_area or area>max_area: continue
        cx,cy = x+w//2, y+h//2
        if any(np.hypot(cx-px, cy-py)<center_dist_thresh for px,py in center_pix):
            continue

        draw.rectangle([(x,y),(x+w,y+h)], outline=outline_color, width=outline_width)

        txt=str(idx); bbox=font.getbbox(txt)
        tw,bh=bbox[2]-bbox[0],bbox[3]-bbox[1]
        draw.text((cx-tw//2, cy-bh//2), txt, fill=text_color, font=font)

        centers.append({
            'index':idx,
            'topleft':(float(round(x / W, 6)), float(round(y / H, 6))),
            'bottomright':(float(round((x+w)/W, 6)),float(round((y+h)/H, 6))),
        })
        center_pix.append((cx,cy)); idx+=1
    return pil_image, centers, center_pix


def get_files_points(las_files_time, resolution):
    logging.info(f"get_files_points 获取点云开始")
    points_time = []
    for fp in las_files_time:
        pts = denoise_las(fp, denose_nb_neighbors, denose_std_ratio)
        points_time.append(pts)
    logging.info(f"get_files_points 获取点云完成")
    return generate_dem_2(points_time, resolution)


def img_resize(img):

    target_height = 2800
    scale = target_height / img.height
    target_width = int(img.width * scale)
    img_resized = img.resize((target_width, target_height), Image.BILINEAR)
    return img_resized

def get_all_files(base_path):

    try:
        las_files_time1 = get_las_files(base_path, 'first_file')
        las_files_time2 = get_las_files(base_path, 'second_file')
        rgb_images_path = base_path + '/tif/rgbimg.tif'
        realland_images_path = base_path + '/tif/realland.tif'
        logging.info(f"时间点1加载了 {len(las_files_time1)} 个LAS文件")
        logging.info(f"时间点2加载了 {len(las_files_time2)} 个LAS文件")
        return las_files_time1, las_files_time2, rgb_images_path, realland_images_path
    except Exception as e:
        logging.info(f"get_all_files 加载文件时出错: {e}")






if  __name__ == '__main__':

    resolution = 0.1           # DEM分辨率，值越小越精细但慢（例如 0.5=每个像素显示0.5平方米的高程值）
    threshold = 0.7            # 高程变化阈值（以米为单位）


    # 设置基础目录
    base_directory = 'E:/A1/A1_0709_0627'  # 根据实际情况修改

    # 加载两个不同时间点的点云文件

    las_files_time1 = get_las_files(base_directory, '0627')
    las_files_time2 = get_las_files(base_directory, '0709')

    logging.info(f"时间点1加载了 {len(las_files_time1)} 个LAS文件")
    logging.info(f"时间点2加载了 {len(las_files_time2)} 个LAS文件")


    # 第一时间点合并后的dem，包括坐标和高程
    gx1, gy1, gz1 = get_files_points(las_files_time1, resolution)
    logging.info(f"gx1.shape={gx1.shape}, gy1.shape={gy1.shape}, gz1.shape={gz1.shape}")

    # 第二时间点合并后的dem，包括坐标和高程
    gx2, gy2, gz2 = get_files_points(las_files_time2, resolution)
    logging.info(f"gx2.shape={gx2.shape}, gy2.shape={gy2.shape}, gz2.shape={gz2.shape}")

    gx3, gy3, gz_diff = dem_diff_overlap(gx1, gy1, gz1, gx2, gy2, gz2)

    logging.info('第一时间DEM', gz1.shape, gx1[0][0], gy1[0][0], gx1[-1][-1], gy1[-1][-1])
    logging.info('第二时间DEM', gz2.shape, gx2[0][0], gy2[0][0], gx2[-1][-1], gy2[-1][-1])
    logging.info('重叠高程差', gz_diff.shape, gx3[0][0], gy3[0][0], gx3[-1][-1], gy3[-1][-1])

    # 只保留超过threshold的高程差（下塌）
    threshold_mask = gz_diff >= threshold

    gz_diff_filtered = np.full_like(gz_diff, np.nan, dtype=float)
    gz_diff_filtered[threshold_mask] = gz_diff[threshold_mask]

    ratio = np.count_nonzero(threshold_mask) / threshold_mask.size
    logging.info(f"过滤后保留超过 {threshold:.2f} m 阈值的点占比：{ratio:.2%}")

    # 过滤掉排队图层中透明部分
    rows_m, cols_m, xs_m, ys_m, dz_m = filter_by_realland_mask(
        gx3, gy3, gz_diff_filtered, f'{base_directory}/realland1.tif'
    )

    # 在底图上叠加下塌点
    blue_img, x_top_left, y_top_left, x_right_bottom, y_right_bottom = overlay_blue_on_rgb(xs_m, ys_m,
                                                                                           f'{base_directory}/rgbimg.tif')

    marked, infos, center_pix = mark_color_blocks_with_number_2(
        blue_img,
        target_rgb=(0, 0, 255), tol=10,
        min_area=500, max_area=15000,
        outline_color=(255, 0, 0), outline_width=2,
        text_color=(0, 255, 0), font_size=30,
        center_dist_thresh=10
    )

    # 保存图像
    marked_pic = os.path.join(base_directory, "final_result_2.png")

    # 泊村指定大小大图片
    image_resize = img_resize(marked)
    image_resize.save(marked_pic)

    logging.info(f"完成红框标记，已保存 {image_resize}")
    logging.info("检测到的框信息：")
    logging.info(f"底图左上角坐标: {x_top_left}, {y_top_left}, 右下角坐标: {x_right_bottom}, {y_right_bottom}")

    # 初始化结果列表
    result = {
        "Rightxy": [  # 可以根据需要补充底图的左上角和右下角坐标
            [x_top_left, y_top_left],
            [x_right_bottom, y_right_bottom]
        ],
        "LatAndLngList": []
    }

    # 计算图像的宽度和高度（地理坐标系）
    image_width = x_right_bottom - x_top_left
    image_height = y_right_bottom - y_top_left

    for item in infos:
        # 从归一化坐标（0-1范围）转换为实际地理坐标
        topleft_lon = x_top_left + item['topleft'][0] * image_width
        topleft_lat = y_top_left + item['topleft'][1] * image_height
        bottomright_lon = x_top_left + item['bottomright'][0] * image_width
        bottomright_lat = y_top_left + item['bottomright'][1] * image_height

        # 计算中心点坐标（左上角和右下角的平均值）
        center_lon = (topleft_lon + bottomright_lon) / 2
        center_lat = (topleft_lat + bottomright_lat) / 2

        # 构造单个坐标对象
        coordinate_obj = {
            "Index": item['index'],  # 序号
            "Center": [center_lon, center_lat],  # 中心点坐标
            "Topleft": [topleft_lon, topleft_lat],  # 左上角坐标
            "Bottomright": [bottomright_lon, bottomright_lat]  # 右下角坐标
        }
        # 添加到结果列表
        result["LatAndLngList"].append(coordinate_obj)
    logging.info(result)





